Confinement control by optical lattices 
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It is shown that the interplay of a confining potential with a periodic potential leads for free 
particles to states spatially confined on a fraction of the total extension of the system. A more 
complex "slicing" of the system can be achieved by increasing the period of the lattice potential. 
These results are especially relevant for fermionic systems, where interaction effects are in general 
strongly reduced for a single species at low temperatures. 
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I. INTRODUCTION 

The study of trapped atomic gases has become a field 
of intense research in the past years. The realization 
of Bose-Einstein condensation (BEC) in trapped dilute 
atomic vapors Ql - Q was the main motivation start- 
ing all the experimental and theoretical research in this 
area. BEC was obtained trapping and evaporatively cool- 
ing bosonic alkali metals. Recently, the possibility of 
trapping and cooling Fermi gases has attracted a lot of 
attention, due to the fact that in the quantum degener- 
acy regime, superfluidity appears within reach How- 
ever, cooling single component Fermi gases up to very 
low temperatures is more difficult than cooling bosonic 
gases since the s-wave collisions arc forbidden for iden- 
tical fermions. On the other hand, single species Fermi 
gases make in this way possible to access experimentally 
an ideal Fermi gas. As shown below, such a simple sys- 
tem can develop a rich behavior by the combination of a 
confining and a lattice potential. 

For atoms confined in a harmonic trap, a case that 
adequately describes most of the experiments realized so 
far (3, a fairly complete theoretical understanding was 
achieved for the one-dimensional (ID) 0, 01 B Hi an d 
two- or three-dimensional |lCt ITU ll2T | single component 
spin polarized trapped Fermi gas, which at very low tem- 
peratures can be considered as a nonintcracting gas. The 
harmonic form of the potential allows obtaining a num- 
ber of exact analytical results for these systems. How- 
ever, these results cannot be extended to incorporate an 
additional lattice potential, a case of increasing interest 
after the experimental realization of a Mott insulator in 
the presence of an optical lattice 01 . A further inter- 
est on the introduction of an optical lattice in fermionic 
systems arises from the possible connections with central 
problems in condensed matter physics [T3 |. 

We analyze here ground state properties of single 
species nonintcracting fermions confined on ID optical 
lattices. These systems are relevant for the understand- 
ing of recent experimental results [TEl fill Il7j , where due 
to the very low temperatures achieved, fermions can be 
considered as nonintcracting particles. On the theoret- 
ical side, the Hamiltonian can be diagonalized numeri- 
cally, which allows to consider any kind of trapping po- 



tential and any number of dimensions for the system. We 
show that the interplay between the lattice and the con- 
fining potential leads in a region of the spectrum to a 
splitting of the system with eigenstates that have a non- 
vanishing weight only in a fraction of the trap. Hence, 
such systems are qualitatively different from the cases 
without the lattice, which have been studied recently 

We also study the noncquilibrium dynamics of the 
fermionic cloud on a lattice. In particular, we study 
the case in which the center of the trap is initially dis- 
placed a small distance. It allows to realize the exis- 
tence of the single particle states confined in a part of 
the trap obtained in the equilibrium case, since for some 
values of the parameters the center-of-mass (cm.) of 
the system oscillates in one side of the trap. With these 
results, we reproduce the experimental observations in 
Rcfs. \13L Ha . fl7| . and complement other theoretical ap- 
proaches to this problem [13, Ha LUJ ■ 

We show that if in addition to the lattice an alternating 
potential is introduced, doubling the original periodicity, 
an additional "slicing" of the system can be achieved. 
The width and number of such regions can be controlled 
in a given energy range by the amplitude of the new mod- 
ulation. By filling these systems with fermions, insulating 
regions may appear, that in the case of an alternating po- 
tential, arc similar to the Mott insulating plateaus of the 
trapped fermionic Hubbard model [2CJ, |21| . In the nonin- 
teracting case it is possible to calculate the local density 
of states, which exhibits the presence of a local gap s in 
the system. In addition, a local compressibility j20t l2lj 
also serves as a local order parameter to characterize the 
insulating regions. This extends the results initially ob- 
tained for the bosonic case |52 , showing that in general, 
the distinction between commensurate and incommensu- 
rate fillings typical in extended solid-state systems is lost 
in the trapped system. 

The presentation is organized as follows: In Sec. II we 
study ID lattices superposed to a confining potential. We 
analyze the generic features valid for any kind of trap- 
ping potential, and focus on fermionic systems. In Sec. 
Ill, an analysis of the noncquilibrium dynamics of the ID 
trapped fermions is presented, and recent experimental 
results reproduced. In Sec. IV we study the case in which 
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an additional alternating potential in introduced, and di- 
cuss analogies and differences with the results obtained 
for the fermionic Hubbard model. In Sec. IV, we extend 
the analysis of Sec. II and IV to two dimensions (2D). 
Finally the conclusions are given in Sec. V. 



II. NONINTERACTING PARTICLES 
CONFINED IN ID OPTICAL LATTICES 

We analyze in this section ID noninteracting systems 
confined by arbitrary potentials when an underlying opti- 
cal lattice is present. We first show results for a harmonic 
confining potential and then discuss the features that are 
generally valid for any other kind of confining potential. 
For definitcness we concentrate on the fermionic case, al- 
though the spectral features are equally valid for bosons, 
since we deal with the noninteracting case. 

In the second quantization language, the Hamiltonian 
describing a confined dilute and ultracold (noninteract- 
ing) gas of single-species fermions, under the influence of 
a ID optical lattice, can be written as 



H 



J dr *t(r) 



2m 



V 2 + V(r) + V (x) 



*(r), (1) 



where ^t(r) anc j \&(r) are the creation and annihilation 
fermionic field operators, respectively The confining po- 
tential is denoted as V(r) = V(x) + V(y) + V(z). We 
analyze in this section the case in which the transversal 
component of the confining potential V(y) + V(z) is very 
strong so that only its lowest energy state is populated, 
and the exited states are not accessible for the given ex- 
perimental setup. Hence, the relevant dynamics of the 
system is restricted to occur in the longitudinal direction 
where the trap is considered to have an arbitrary power 
a, V(x) = V a x a . In Eq. Q, V (x) = V cos 2 (kx) de- 
scribes the potential generated by a ID optical lattice. 
The wave vector k = 2ir/\ is determined by the wave- 
length A of the laser beam. (The lattice spacing is then 
a = A/2.) Assuming the atoms to be at the lowest vi- 
brational level in each site, the fermionic field operators 
can be expanded in single band Wannicr functions 4>i(x), 
= ^Cj^i(a;), and from Eq. one obtains the 
single band Hamiltonian 



H 



eje, 



i+1 



H.c 



X; m 



(2) 



where c\ and c, are creation and annihilation operators, 
respectively, for a spin polarized fermion on site i, the lo- 
cal density is n, = cjc i; and xi measures the positions of 
the sites in the trap (xi = ia with —N/2 + 1 < i < N/2, 
N being the number of lattice sites). The hopping 
parameter is denoted by t, which for Vq S> E r can 
be written in terms of the experimental parameters as 
t = A/y^(V /E r ) 3 ^E r e- 2 V v o/ E - [53, where the recoil 
energy of the atoms (with mass m) is E r = H 2 k 2 /2m. 



The total number of spin polarized fermions in the sys- 
tem is denoted by Nf. We diagonalize the Hamiltonian 
numerically, and consider the cases in which all particles 
are confined. 




2000 4000 6000 8000 10000 
level number 

FIG. 1: Single particle spectrum (a) and level spacing (b) for 
a system with N = 10000 and Via 2 = 3 x 10~ 7 £. Energies are 
measured in units of t. For the explanation of the arrows see 
text. 



Results obtained for the single particle spectrum of a 
system confined by a harmonic potential are presented in 
Fig. [If a). The spectrum is clearly different from the usual 
straight line in the absence of a lattice. It is possible to 
see that in Fig. ^a) the spectrum can be divided in two 
regions according to the behavior of the energy as a func- 
tion of the level number. An arrow is introduced where a 
change in the curvature is observed. More detailed infor- 
mation can be obtained by considering the level spacing 
as a function of the level number [Fig. IHb)]. There it 
can be seen that in the low energy part of the spectrum 
(region A), the level spacing decreases slowly with in- 
creasing level number, in contrast to the case without the 
lattice in which the level spacing is constant. However, 
at the point signaled with the first arrow, a qualitative 
change in the single particle spectrum occurs, character- 
ized by an oscillating behavior of the level spacing. The 
part with values of the level spacing increasing with the 
level number corresponds to odd level numbers and the 
one with a level spacing that decreases up to zero cor- 
responds to even level numbers. That is, a degeneracy 
sets in that continues up to the point signaled with the 
second arrow, where a new change in the behavior of the 
level spacing shows up. The region beyond the second 
arrow corresponds to deconfined states, which are of no 
interest since experimentally they are associated to par- 
ticles that scape from the trap (which in the system of 
Fig. U] has 10000 lattice sites). 
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In the lowest part of the spectrum of Hamiltonian J2j , 
the eigenfunctions are essentially the harmonic oscilla- 
tor (HO) orbitals in the absence of a lattice. This is 
shown in Fig. [2{a) f° r the first and the second eigen- 
functions of Eq. and the same parameters of Fig. 
n These orbitals are perfectly scalable independently 
of the size of the system and of the ratio between V2 
and t. It is only needed to consider that the usual HO 
characteristic length R ~ (mw) (without the lat- 
tice) is given in terms of the lattice parameters through 

R ~ (V2/ta 2 ) with the effective mass m ~ (to 2 ) 1 
for very low energies, so that the scaled orbitals are given 

1 /2 

by <p = (R/a) 4> where <j> are the HO orbitals with 
the lattice, i.e., the same relation as for the HO with- 
out the lattice holds for the lowest energy orbitals with 
the periodic potential. This implies that very dilute sys- 
tems behave similarly to continuous systems, which have 
been already discussed in the literature so that we do 
not present any further analysis on them. The Nf os- 
cillations in density profiles and momentum distribution 
function (MDF), and other mentioned characteristics of 
the ID trapped system without the lattice 0,13 are easily 
obtained in this case. 
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FIG. 2: Scaled HO orbitals in the presence of a lattice for N = 
10000 and V 2 a 2 = 3 x 10 -7 t. (a) First (continuous line) and 
second (dashed line) HO orbitals, (b) HO orbitals 5000 (only 
different from zero for negative x) and 5001 (only different 
from zero for positive a;). In (a) the positions are given in 

units of the HO length R = (V2/ia 2 ) (for an explanation, 
see text), and in (b) in units of the lattice constant a. 

A qualitative difference between the cases of the trap 
with and without a lattice starts for levels in region B. 
Once the degeneracy appears in the spectrum, the cor- 
responding eigenfunctions of the degenerate levels start 
having zero weight in the middle of the trap, and for 
higher levels the regions over which the weight is zero in- 
creases. As an example, we show in Fig.[2Jb) two normal- 
ized eigenfunctions belonging to region B in Fig. ^ The 
cases depicted correspond to the normalized eigenfunc- 
tions 5000 (that is only different from zero for negative 
values of x) and 5001 (only different from zero for posi- 
tive x), for the same parameters of Fig.n( m principle a 
lineal combination of these two eigenfunctions could have 
been the solution since the level is degenerated). Hence, 



particles in these states are confined to a fraction of the 
trap, showing that the combination of both a confining 
and a periodic potential lead to features not present ei- 
ther in the purely confined case without a lattice or in the 
case of a purely periodic potential. Furthermore, since we 
arc dealing with a nonintcracting case, such features are 
common to both fermions and bosons. However, in the 
case of fermions, it is easy to understand the reason for 
such effects, as discussed next. 

Figure |3f a) shows density profiles of fermions when 
the number of particles in the trap is increased. In one 
case (Nf = 4500) the Fermi energy lies just below the 
level marked with an arrow in Fig. ^ A second curve 
(Nf = 4651) corresponds to the case where the central 
site reaches a density n = 1, and in the other case (Nf = 
5001), the Fermi energy lies at the value corresponding 
to the levels depicted in Fig.[5Jb). The positions in the 
trap are normalized in terms of the characteristic length 
for a trapped system when a lattice is present, which is 
given by EE El 

C = (v 2 /ty 1/2 . (3) 

When the Fermi energy approaches the level where de- 
generacy sets in, the density of the system approaches 
n — 1 in the middle of the trap, and at the filling point 
where the degeneracy appears in the spectrum, the den- 
sity in the middle of the trap is equal to one, so that 
an insulating region appears in the middle of the system. 
Increasing the filling of the system increases the region 
over which this insulator extends. Hence, due to Pauli 
principle, the eigenfunctions of such levels cannot extend 
over the insulating region, and for the same reason, the 
region over which the weight is zero increases for higher 
levels. The local insulator with n = 1 has zero variance 
of the density and, it is incompressible, a property that 
could be tested experimentally by using a local probe. 
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FIG. 3: Density profiles (a) and normalized MDF (b) for 
Nf — 4500 (thick continuous line), Nf — 4651 (dashed line), 
and Nf = 5001 (think continuous line) for a system with 
N = 10000 and V 2 a 2 = 3 x 10~ 7 t. In (a) the positions are 
given in units of the characteristic length £, and in (b) the 
momentum is normalized by the lattice constant a. 

However, since the system we are considering is a 
nonintcracting one, the confined states discussed above 
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should be also present in the case of bosons, where the ar- 
gument about the filling would not be valid anymore. It is 
therefore desirable to also understand such features from 
a single particle perspective [24J. We first notice that 
the point at which degeneracy appears is at an energy 4i 
above the lowest level (Eq) [see Fig. ^a)], correspond- 
ing to the bandwidth for the periodic potential. Such an 
energy is reached when the Bragg condition is fulfilled 
and in the case of the tight-binding system we are con- 
sidering, when all the available states are exhausted. Let 
us next consider the case depicted in Fig. |2fb). There, 
the energy level corresponding to the wave functions is 
-S5001 — Eq = 4.2176i, that is to a good approximation 
At + V 2 x\ for V 2 a 2 = 3 x l(T 7 i and xi = 697a the 
inner point where the wave functions drop to a value 
~ 1CP 5 . Therefore, the inner turning point corresponds 
to the Bragg condition, whereas for the outer turning 
point (x2 = 3770a, again for the same drop of the wave 
function), we have that E 500 i — E ~ V2X2, i.e., the clas- 
sical turning point corresponding to the harmonic po- 
tential, as expected for such a high level. Hence, B rag g 
scattering as in the well known Bloch oscillations |25j . 
and the trapping potential combine to produce the con- 
finement discussed here. 

Further confirmation of the argument above can be ob- 
tained by considering the MDF, a quantity also accessible 
in time of flight experiments [l3j. Due to the presence 
of a lattice, it is a periodic function in the reciprocal lat- 
tice |26j and it is symmetric with respect to k = 0, so 
that we study it in the first Brillouin zone in the region 
[0, ir /a]. In addition we normalize the MDF to unity 
at k = (/ife=o = 1). For the fermionic case, it can be 
seen that it always has a region with nk = if the in- 
sulating phase is not present in the trap, and this region 
disappears as soon as the insulator appears in the mid- 
dle of the system. More precisely, Fig. 02b) shows that 
at the filling when the site in the middle reaches n = 1, 
also the momentum k = ir/a is reached, such that the 
Bragg condition is fulfilled for the first time, confirming 
the discussion above. When further sites reach a density 
n = 1, n^—^/a increases accordingly. Then the formation 
of the local insulator in the system can be tested ex- 
perimentally observing the occupation of the states with 
momenta k = ±ir/a. 

Since for different systems sizes and number of parti- 
cles, potentials with different curvatures have to be con- 
sidered, it is important to determine the filling at 
which the insulator appears in the middle of the trap 
as a function of the curvature of the harmonic confining 
potential. This question was already answered for the 
interacting case (Hubbard model) in Refs. where 
we determined the phase diagram. There we showed that 
if a dimcnsionless characteristic density p is defined as 
p = Nfd/C, then its value when the insulating regions 
(Mott insulating and band insulating in the interacting 
case) appear in the system is always constant for any 
value of Vz/t at a given value of U/t (within error bars 
there), so that Nf ~ (/a. However, in Refs. [H ^ 



we were able to check this only up to 150 lattice sites 
and fillings up to the same order, whereas here we ex- 
tend those results to much larger systems. In Fig. Ufa) 
we show in a log- log scale how Nj depends on V2 /t over 
three decades on the total filling. In our fit the slope 
of the curve is —0.500 (with 0.04 percent of error), as 
expected on the basis of Eq. © . The critical character- 
istic density pc = N^a (V^/i) 1 ^ 2 at which the insulating 

region appears is pc = e^, with j3 = 0.986 (with 0.3 per- 
cent of error), which is curiously rather close to the basis 
of the natural logarithms. 
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FIG. 4: Total filling in the trap needed for the formation of 
the insulator as a function of the curvature of the confining 
potential, (a) Harmonic potential, (b) Potential with a power 
a = 6. 

For systems with other powers for the confining po- 
tential it is only needed to define the appropriate dimcn- 
sionless characteristic density p — Nfa (Va/t) 1 ^ 01 , and 
determine its value at the point where the insulator ap- 
pears. In Fig. Efb) wc show in another log-log plot how 
Nj? depends on the curvature of a confining potential 
with power six (Ve/t). As anticipated, we obtain that 
the slope of the curve is 1/6 (with 0.01 percent of error) 
in this case and the characteristic density for the forma- 
tion of the insulator is pc = 2.09. Fina lly, we should 
mention that it was already shown in Ref. |2l| that keep- 
ing constant the characteristic density but changing the 
curvature of the confining potential and the total fill- 
ing in the trap, the density profiles as a function of the 
normalized coordinate and the normalized MDF remain 
unchanged. 

In general, for arbitrary confining potentials the same 
features discussed previously for the harmonic case are 
valid. The spectrum and level spacing behave in a dif- 
ferent way depending on the power of the confining po- 
tential, but always at a certain level number degeneracy 
appears in the single particle spectrum and it corresponds 
to the formation of an insulator in the middle of the sys- 
tem for the corresponding filling. In Fig. we show the 
single particle spectrum [Fig. Efa)] and the correspond- 
ing level spacing [Fig.[5{b)] for a confining potential with 
power a = 6, where the features mentioned previously 
are evident. The arrow in the inset of Fig. shows the 
level at which degeneracy sets in, very much in the same 
way as in the harmonic case. 
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7000 14000 21000 28000 

level number 

FIG. 5: Single particle spectrum (a) and level spacing (b) 
for a system with a confining potential with power a = 6, 
N = 28000, and V 6 a 6 = 3 x 10 7 t. Energies are measured 
in units of t. Expanded view in (b) shows the first part of 
the level spacing. An arrow was introduced in the inset for 
signaling the level at which the degeneracy appears in the 
spectrum. 

We close this section by considering the pair distribu- 
tion function. This quantity not only reflects the conse- 
quences of Pauli's exclusion principle but clearly charac- 
terizes the insulating region. In the presence of a lattice 
the pair distribution function can be written as 

Pi i = (n i )(n j )-ti jl (4) 

where pij = (cjcj) is the fermionic one-particle density 
matrix. 

In Fig. we show as intensity plots the pair distribu- 
tion function for systems with TV = 1000 lattice sites and 
V^a 2 = 3 x 10 -5 t. Figure Ufa) corresponds to the case 
with Nf = 300 fcrmions, where the systems is completely 
metallic, whereas Fig. E[b) corresponds to Nf = 600 
fcrmions, an insulating region appears in the middle of 
the trap. Apart from the depression along the diagonal 
that reveals the consequences of Pauli's exclusion princi- 
ple, a clear distinction between the purely metallic case 
and the one with an insulating region can be seen. In- 
side the insulating region, the density matrix becomes 
diagonal, such that Pij = 1 for i ^ j and Pa =0. 

III. OSCILLATIONS OF FERMIONS IN A ID 
LATTICE 

Recent experiments have realized single species non- 
interacting fermions in ID optical lattices 0, ITfl [r^| . 




FIG. 6: (Color online) Intensity plots of the pair distribution 
function as a function of the normalized position for (a) Nf = 
300 (only a metal is present in the trap) and (b) Nf = 600 
(a insulating region is formed in the middle of the trap). The 
number of lattice sites is N = 1000 and V 2 a 2 = 3 x 10~ 5 t. 



Transport studies in such systems revealed that under 
certain conditions a sudden displacement of the trap cen- 
ter is followed by oscillations of the cm. of the fermionic 
cloud in one side of the trap. This is in contrast to the 
system without the lattice where the cm. oscillates, as 
expected, around the potential minimum 

HEmm Al- 
though the experimental system is not a true ID system, 
due to the strong transversal confinement the relevant 
motion of the particles occurs in the longitudinal direc- 
tion. Hence, in order to qualitatively understand the ob- 
served behavior one can analyze the ideal ID case. Given 
the results discussed in the previous section one expects 
the displaced oscillation of the cm. to appear when, 
due to the initial displacement of the trap, particles that 
where located in region A of the spectrum in Fig. ^b) 
are moved into region B so that Bragg conditions are ful- 
filled. Then the particles get trapped in one side of the 
system [Fig.^b)]. 

Figure [7] shows exact results obtained for the cm. dy- 
namics of 1000 fcrmions in a trap with N = 3000 when 
its center is suddenly displaced 200 lattice sites, r de- 
notes the real time variable. The relation between the 
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confining potential (Va) and the hopping parameter (t) 
is increased in order to fulfill the Bragg conditions. This 
is equivalent in experiments to increase the curvature of 
the confining potential keeping constant the depth of the 
lattice, which leads to an increase of the frequency of 
the oscillation as shown in Fig. It is also equivalent 
to increase the depth of the lattice keeping the confining 
potential constant, but then our plots in Fig.[7]should be 
interpreted with care since there we normalize the time 
variable by the hopping parameter, which changes in the 
latter case. In Fig. [7\ (dashed line) we show results for 
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FIG. 7: Evolution of the cm. (x C m) of 1000 confined fermions 
when the center of the trap is suddenly displaced 200 lattice 
sites (x ), for V^a 1 = 2 x 10~ o6 i (dashed line), and V^a? = 
6 x W~ 06 t (continuous line). 

the case where the cm. of the cloud oscillates around the 
minimum of energy of the trap since no Bragg conditions 
are fulfilled. This can be seen in the MDF [Fig. Eta)] 
where at any time no particles have k = ±7r/a. Figure [7| 
also shows that a damping of the oscillation of the cm. 
occurs. This is due to the nontrivial dispersion relation 
in a lattice tu = —2tcoska, which makes the frequency 
of oscillation of the particles dependent on their energies, 
leading to dephasing. In order to reduce the damping, 
fermions should populate after the initial displacement 
only levels with energies close to the bottom of the band 
in a lattice, so that the quadratic approximation is valid 
for gfe. (Notice that this is not generally fullfiled even if 
the initial displacement is small.) 

Increasing the relation V^/t makes that some particles 
start to fulfill the Bragg conditions so that the center of 
oscillations of the cloud depart from the middle of the 
trap. In Fig. [7| (continuous line) we show a case where 
the cm. never crosses the center of the trap. The MDF 
corresponding to this case, at three different times, is 
displayed in Fig. IBtb). There it can be seen that ini- 
tially (t = 0) no Bragg conditions are satisfied in the 
system, and that some time after the initial displacement 
the Bragg conditions are fulfilled (r = 250h/t). Finally, 
we also show the MDF long time after the initial displace- 
ment of the trap (r = 12500ft/t), when the oscillations 
of the cm. are completely damped and the MDF is ap- 
proximately symmetric around k = 0. 




FIG. 8: MDF of 1000 trapped fermions at three different times 
after displacing the trap 200 lattice sites, for V2C1 2 = 2x 10 _06 t 
(a), and V 2 a 2 = 6 x 10~ 06 t (b). The times are r = (dotted 
line), r = 250ft/t (dashed line), and r = 12500ft/t (continuous 
line) . 



IV. DOUBLING THE PERIODICITY 

In this section we study the consequences of enlarging 
the periodicity in the lattice. For this purpose we in- 
troduce an alternating potential, and the Hamiltonian of 
the system can be written as 



H = -*E( c i c *+i +H - c -) + ^E^) c 

i i 



(5) 



where the last term represents the oscillating potential 
and V a its strength. The purpose of introducing an al- 
ternating potential in the trapped system is twofold. For 
fermionic systems, the increase of the numbers of sites 
per unit cell leads to the possibility of creating insulat- 
ing states (band insulators in the unconfincd case) for 
commensurate fillings. On the other hand, by changing 
the periodicity, new Bragg conditions are introduced, giv- 
ing the possibility of further control on the confinement 
discussed in the previous sections. 

Figure [5] shows how the density profiles evolve in a 
harmonic trap when the total filling is increased. Since 
the density oscillates due to the alternating potential, we 
made two different plots for the odd [negative value of 
the alternating potential, Fig. Eta)] and even [positive 
value of the alternating potential, Fig. Elb)] sites. Each 
of the plots in Fig.[3]is very similar to the evolution of the 
density profiles already shown for the trapped Hubbard 
model 123. The only difference is that in Figs. [Ja) 
andOJb) the plateaus with n ^ 1 have densities different 
between themselves and different from n = 0.5, which 
would be the density of one component of the spin polar- 
ized fermions in the Mott insulating phase of the Hub- 
bard model. In the flat regions of Fig. EJ both even and 
odd sites have the same densities than the corresponding 
sites in the periodic case at half filling for the same value 
of the alternating potential, so that it is expected that 
they correspond to local insulating phases. 
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FIG. 9: Evolution of the local density in a harmonic trap as a function of the position and increasing total number of fermions 
when an alternating potential V a = 0.5t is present, (a) Odd sites, (b) Even sites. The system has 1000 lattice sites and 
V 2 a 2 = 3 x 10~ 5 t. 



In Figs. II Uf a) and HUf b') we show the single particle 
spectrum and the level spacing respectively for the same 
parameters of Fig. |5J Although in this case the level 
spacing exhibits a more complicated structure, an imme- 
diate identification between the regions signaled in Fig. 
IToTh) between arrows and different fillings in Fig. can 
be done. (A) corresponds to the fillings in Fig. where 
only a metallic phase appears in the trap, (B) to the fill- 
ings where the first plateau is present in Fig. (C) to 
the fillings where a metallic phase develops in the middle 
of the trap and it is surrounded by insulating regions, 
and (D) to the fillings in Fig. where the insulator with 
n = 1 appears in the center of the system. The region 
after the last arrow in Fig. llUf b) corresponds to decon- 
fined states. Notice that the level spacing in regions (C) 
and (D) shows a behavior that was not present in Fig.^ 

In order to understand the complex behavior of the 
level spacing we study, as in the previous section, the 
eigenfunctions of the system shown in Fig.^| The eigen- 
functions corresponding to region A in Fig. ^b) behave 
as expected for a metallic phase, where the combina- 
tion of the alternating and confining potentials generates 
a different modulation than the one studied in Sec. II, 
but without qualitative differences. In the second region 
of the spectrum (region B) the eigenfunctions have zero 
weight in the middle of the trap, exactly like in the insu- 
lator discussed in Sec. II. In region C there is, as pointed 
out above, a new feature since in this case it is possible 
to obtain a metallic region surrounded by an insulating 
one. This is reflected by the eigenfunctions shown in Fig. 
Illf a). where one of the eigenfunctions is nonzero only 
inside the local insulating phase (continuous line), and 
the other is nonzero only outside the insulating phase 
(dashed line) , the energy levels associated with the latter 
ones are degenerated. For the region D the situation is 
similar but in this case the system is divided in four parts 
because of the existence of the insulator with n = 1 in 
the middle of the trap and the insulator between the two 
metallic phases. This implies that all the levels are de- 




200 400 600 800 1000 

level number 



FIG. 10: Single particle spectrum (a) and level spacing (b) 
for a system with an alternating potential V a = 0.5t and with 
N = 1000 and V 2 a 2 = 3 x 10~ 5 i. Energies are measured in 
units of t. For the explanation of the regions between the 
arrows, see text. 



generate in region D, and the particles are located cither 
between both insulating regions or outside the outermost 
one, as shown in Fig. Illf b). 

As in the previous section, the spectral features dis- 
cussed here are equally valid for fermions as well as for 
bosons. Up to now we discussed the "slicing" of the sys- 
tems only in terms of fermions and based on the appear- 
ance of insulating regions along the system. As before, it 
would be also here desirable to understand the appear- 
ance of forbidden regions in space in terms of a single 
particle picture. We show now that with the introduction 
of new Bragg conditions, due to the altered periodicity, 
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FIG. 11: Eigenfunctions for a trapped system with an al- 
ternating potential (V a = 0.5t) for N = 1000 and V 2 a 2 = 
3 x 10~ 5 i. The eigenfunctions correspond to the levels: 269 
(a) (continuous line), 270 (a) (dashed line), 501 (b) (continu- 
ous line), and 502 (b) dashed line. 

the "slicing" of the system can be explained in an anal- 
ogous way as in the previous section. In the unconfincd 
case, the doubling of the periodicity creates new Bragg 
conditions at k = ±7r/2o, such that an energy gap 214 
appears. Figure llUf a) shows that in the confined case 
the spectrum is continuous (in the sense that the level 
spacing is much smaller than 214), so that the imprint 
of the gap can be seen only in the local density of states 

Ni(u) = -TmGu(u) , (6) 

where Gij(u) is the one-particle Green's function |27j . 
which in this case can be easily computed. 
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FIG. 12: Density profile per unit cell (now containing two 
contiguous lattice sites) for a trap with N = 10000, Via 2 = 
3 x 10~ 7 t, and Nj = 2000. The insets show the density of 
states per unit cell for two points in the profile. The arrows 
in the insets signal the Fermi energy for the selected filling. 

The insets in Fig. ^| show the density of states per 
unit cell (now containing two lattice points) for two dif- 
ferent positions along the density profile. The downward 
arrows in each inset corresponds to the location of the 
Fermi energy. The inset at the left corresponds to a sit- 
uation where the Fermi energy goes through the lowest 
band, whereas the inset at the right belongs to sites in 



the middle of an insulating region. As expected, in this 
latter case, the Fermi energy lies inside the gap. The size 
of the gap is to a high degree of accuracy 214 for the site 
in the middle of the trap, but slightly less on the sides. 
Therefore, again the same arguments as before can be 
used, but inste ad of 4i, the width for each band is given 
by y/ At 2 + 14 2 — 14 — 1.56t in our case. Without repeat- 
ing the detailed discussion in the previous section, we 
can understand the confinement in Fig. Illf a) as follows. 
Level 269 has an energy that for sites in the middle of the 
trap falls in the middle of the upper band, while for level 
270 (the same value of energy), passes through the low- 
est band. In fact, the density of states shown in Fig. 1121 
can be viewed as approximately shifted by V2X 2 , count- 
ing the sites from the middle. Finally, levels in Fig. lllf b) 
correspond to the case where in the middle of the trap 
they fall beyond the highest band, then going outwards, 
they fall in the middle of the highest band, and further 
outside, they fall in the middle of the lowest band. 




-3 -2 -1 1 2 3 1/4 nQ 3H/4 JI 

tlx, ka 



FIG. 13: Density profiles (left) and their normalized MDF 
(right) for N f = 150 (a),(b), 200 (c),(d), 350 (e),(f), 600 
(g),(h) and N = 1000, V a = 0.5t, V 2 a 2 = 3 x 10~ 5 t. 

Again, as in Sec. II, one can follow the same reasoning 
by considering the MDF. Due to the new periodicity it 
displays new features, associated with the fact that in- 
creasing the periodicity the Brillouin zone is decreased, 
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and in the present case a second Brillouin zone is visi- 
ble. In Fig. El we show the density profiles (left) that 
characterize the four different situations present in Fig. 
They correspond to fillings of the trap in the four re- 
gions of the single particle spectrum discussed previously 
in Fig. El Notice that in the figures we included all the 
odd and even points in the density profiles. We plotted 
as horizontal dashed lines the values of the densities in 
the band insulating phase of the periodic system for the 
odd and even sites (from top to bottom respectively), so 
that it can be seen where are located the local insulating 
phases in the trap. The corresponding normalized MDF 
are presented in Fig. El (right). 

In Figs. ll3f a) and llflf b) it is possible to see that when 
only the metallic phase is present in the trap, in the MDF 
an additional structure appears after 7r/2, correspond- 
ing to the contribution from the second Brillouin zone. 
When a first insulating phase is reached, by coming to 
the top of the lowest band, k = n/2a is reached, and in- 
creasing the fillings of the system beyond that point, the 
dip around k = it /2a disappears [Figs. ITST c'l and !13f d)]. 
On adding more particles to the system a metallic phase 
appears inside the insulating plateau, [Figs. I13f e'l and 
I13f f1]. When this metallic phase widens, decreasing the 
size of the insulating phase, nu starts to be similar to the 
rik of the pure metallic phase in the system without the 
alternating potential. Increasing even further the filling 
of the system, when the trivial insulator (n = 1) appears 
in the center of the trap, the tail with very small values 
of rife disappears (like in the system without the alter- 
nating potential the region with zero also disappears 
[Fig. Efb)] and the further increase of the filling in the 
system makes nj, flatter [Figs. I'Hfl' h') and I13f i)] . 

Up to this point, several quantities, like density pro- 
file, pair distribution function, or local density of states 
were taken as evidence for the existence of an insulating 
phase, but a quantitative criterion in the sense of an order 
parameter to characterize the phases was not given. As 
shown already in the case of the Hubbard model [2(], El , 
it is possible to define a local compressibility: 

K i = XM+j > ( 7 ) 

\j\<HVa) 

where 

Xi,j = (mnj) - (rii) (rij) (8) 

is the density-density correlation function, and £(V a ) ~ 
b£(V a ), with £,(V a ) the correlation length of m the pe- 
riodic system at half-filling for the given value of V a . As a 
consequence of the band gap opened in the band insulat- 
ing phase at half filling in the periodic system, density- 
density correlations decay exponentially and there £(V a ) 
can be determined. The parameter b is considered b ~ 10 
(see discussion in Ref. |21|). When this definition is ap- 
plied to the different fillings of Fig. the local compress- 
ibility is zero in the insulating regions and nonzero in 
the metallic phases. The local quantum critical behavior 



found in Ref. |20j at the transition between the metallic 
and Mott insulating phase is not present here since there 
are no interactions between the particles that could gen- 
erate quantum criticality. 

Finally we analyze the phase diagram for these sys- 
tems. It can be generically described by the characteristic 
density p, like the Hubbard model and the noninteracting 
case in Sec. II. In Fig. 1141 we show two phase diagrams 
for two different values of the curvature of the confining 
potential, V 2 a 2 = 3 x 10~H and V 2 a 2 = 3x 10 _4 i. There 
it can be seen that although there is one order of magni- 
tude between the curvatures of the confining potentials, 
the phase diagrams are one on top of the other, the small 
differences are only due to the finite number of particles 
which make the changes in p discrete. Therefore, the 
characteristic density allows us to compare systems with 
different curvatures of the confining potential, number of 
particles and sizes. In addition we checked that keeping 
the characteristic density constant for a given value of V a , 
the density profiles as a function of the normalized coor- 
dinates and the normalized MDF do not change when the 
number of particles or the curvature of the confining po- 
tential are changed in the system, as we already pointed 
out for the case without alternating potential (Sec. II). 




2 4 6 8 10 

Va/t 

FIG. 14: Phase diagram for systems with V20 2 = 3 x 10 _5 t 
(V) and V 2 a 2 = 3 x 10~H (CO- The different phases are 
explained in the text. 

The different phases present in Fig. El arc (A) a pure 
metallic phase, (B) an insulator in the middle of the trap 
surrounded by a metallic phase, (C) a metallic intrusion 
in the middle of the insulator, (D) an insulator with n — 1 
in the center of the trap surrounded by a metal, an insu- 
lator and the always present external metallic phase. For 
very small values of the alternating potential (V a ), phase 
B is not present, and the insulator surrounding the metal- 
lic phase in the center of the trap disappears leaving a 
full metallic phase at the very beginning of phase C. Sim- 
ilarly, the insulator with n = 1 is surrounded only by a 
metallic phase (at the very beginning of phase D). How- 
ever, these regions are very small in the phase diagram 
and we did not include them. At V a = the results of 
Sec. II are recovered since up to p = 2.68 the system is a 
pure metal and for higher characteristic densities there is 
an insulator with n = 1 surrounded by metallic phases. 
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There is one important difference between the phase 
diagram in Fig. 1141 and the one of the trapped Hubbard 
model [2(1 l2ll | . In Fig. El the boundary between regions 
A and B changes appreciably when the value on the alter- 
nating potential is increased while in the Hubbard model 
case (for the values of U that we simulated) it was found 
independent of the value of U. This is possibly due to 
the fact that for the alternating potential, increasing V a 
changes the local densities of the insulating phase while 
in the local Mott insulating phase the density is always 
constant independently of the value of U. 



V. THE 2D SYSTEM 

In this section we extend to 2D the results obtained in 
previous sections for the ID case. The Hamiltonian in 
this case can be written as 

H = -* E ( c k- + H - c -)+E ( v <**,* x ™* + v <*y,yy? y ) n i> 

(9) 

where (a?i,j/i) are the coordinates of the site i, and 
refers to nearest neighbors. The last term in Eq. @ 
allows to consider different strengths V aa ,xi Va .y and 
powers a x , a y of the confining potential in the x, y di- 



rections. We call in what follows N x and N y the number 
of lattice sites in the x and y directions, respectively. 

In Fig. 1151 we show the single particle spectrum and 
its corresponding level spacing for a system with N x = 
N y = 100 lattice sites confined by a harmonic potential 
with V2. x a 2 = V2. y a 2 = 5 x 10~ 3 i. Figure IH)T b) shows 
that degeneracy sets in at the very beginning of the ex- 



pectrum, and this is because of the symmetries of the 
square lattice. In 2D the formation of the insulator in 
the middle of the trap does not generate additional de- 
generacies in the system since it does not split the trap 
in independent identical parts. Then in contrast to the 
ID case no information of its formation can be obtained 
from the level spacing. 
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FIG. 15: Single particle spectrum (a) and level spacing (b) 



for a system with N x 



Ny 



100, and V^cT = V 2 , y a z 



5 x 10 t. Energies are measured in units of t. 




FIG. 16: (Color online) 2D density profiles (a),(d), 2D normalized MDF (b),(e), and intensity plots of the 2D normalized MDF 
(c),(f), for N f = 1000 (a)-(c) and N f = 4000 (d)-(f) fermions in a system with N x = N v = 100 and V 2 , x a 2 = V 2 , y a 2 = 5 x W~ 3 t. 
The color scale in (c) and (f ) is the same as in Fig. |S] 
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Two density profiles and their corresponding normal- 
ized MDF for Nf = 1000 and N f = 4000, and the 
same trap parameters of Fig. El ar e presented in Fig. 
1161 The x and y coordinates in the trap are normal- 
ized by the characteristic lengths ( x = (V^/t) and 

Cy = (V2, y /t) , respectively. The density profile in 
Fig. I16f a) corresponds to a pure metallic phase in the 
2D trap, the corresponding MDF [Fig. Ilfif b)] is smooth 
and for some momenta it is possible to see that n k = 
like in the ID case. When the filling of the system is 
increased, the insulator appears in the middle of the trap 
[Fig. nUT d - )] and all the regions with n k = 0, present in 
the pure metallic phase, disappear from the MDF [Fig. 
llGf el]. Figs. I16f c) and ll6f f) show as intensity plots the 
normalized MDF of Figs. EJb) and llfit e). 

In 2D it is possible to define a dimcnsionless charac- 
teristic density as p = Nfa 2 /( x ( y . Also in this case it 
has always the same value when the insulator appears 
in the middle of the system, independently of the values 
and relations between V^ x and V2 y - The density profiles 
as function of the normalized coordinates and the MDF 
remain unchanged when the characteristic density is kept 
constant and the values and relations between V% x and 
Viy are changed (in the thermodynamic limit they have 
the same form shown in Fig. I16|) . This implies that the 
results shown in Fig. 1161 for a symmetric trap do not 
change for an asymmetric trap with the same character- 
istic density. The value of the characteristic density for 
the formation of the insulator in a harmonic 2D trap is 
pc ~ 13.5. 

The addition of the alternating potential leads to re- 
sults similar to those presented in the ID case. Four 
density profiles showing the possible local phases in the 
2D trap, and intensity plots of their corresponding MDF 
are shown in Fig. El 111 the pure metallic case [Figs. 
I17f a) and ll7f b1] the additional structure in the MDF for 
k x ,k y > ir/2a, due to the increase of the periodicity, is 
present. This structure also disappears when the insula- 
tor appears in the middle of the system [Figs. I17f c') and 
I17f d)]. Increasing the filling a new metallic phase ap- 
pears m the center of the trap [Figs. Etc) andfTTTd)]. 
For the highest filling the insulator with n = 1 develops 
in the middle of the trap and the MDF becomes flatter 
with nfc everywhere. In a way similar to the ID case 
there are confined states in the radial direction, i.e., par- 
ticles are confined in rings around the center of the trap, 
and they can be explained in terms of Bragg conditions. 
The phase diagram of the 2D case is also similar to the 
one in the ID case and is not discussed here. 

We close this section by considering the local compress- 
ibility. As it was mentioned in the analysis of the ID case 
with the alternating potential, this quantity is zero in the 
insulating phases (like in the Mott insulating phases of 
the trapped Hubbard model). In the 2D case we extend 
the definition given by Eq. to 

K i^iy = XixiySx+jx'iy+jy 1 (10) 

\j*,3v\<l(Vo.) 




FIG. 17: (Color online) Density per unit cell (a),(c),(e),(g) 
and intensity plots of the normalized MDF (b),(d),(f),(h) pro- 
files for N f = 200 (a),(b), 800 (c),(d), 1100 (e),(f), and 3000 
(h),(i) in a system with N x = N y = 100, V^a 2 = V2, y a 2 = 
5 x 10 -3 t, and V a = t. The color scale in the intensity plots 
of the normalized MDF is the same as in Fig. |H] 



where 

Xi.iyj.jy = ("Wv^-Jv) - ( n «*0 ( n Mv) ' (11) 

is the density-density correlation function in 2D and 
i(V a ) — b£(V a ). In this case it is also possible to de- 
termine £(Va) in the insulating phase of the 2D periodic 
case (at half filling) and apply the new definition to the 
2D trap (as for the ID case b ~ 10). The results ob- 
tained for the same parameters of Fig. El are presented 
in Fig. El There it can be seen that the rings of local 
insulators in Fig. 1171 are represented in Fig. 1181 by rings 
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FIG. 18: (Color online) Intensity plots of the local compress- 
ibility as a function of the normalized coordinates for systems 



with N x = N y = 100, 



Vi, 



= V 2 , 



= 5 x lQ~ s t and V a = t. 



The fillings are the same as that in Fig. El (a) Nf = 200, (b) 
N f = 800, (c) Nf = 800, and (d) N f = 1100. Black color 
means zero compressibility. 



of incompressible regions (black rings) so that this defi- 
nition works also perfectly in the 2D case, and the local 
compressibility should be a relevant quantity to charac- 
terize the local Mott insulating phases also in the trapped 
2D Hubbard model. 



VI. CONCLUSIONS 

We performed a detailed analysis of nonintcracting sys- 
tems focusing on the consequences of the combination of 
a confining and a periodic potential. It leads to a con- 
finement of particles in a fraction of the available system 
size. This confinement is directly related to the forma- 



tion of insulating regions in the case of fermionic systems. 
Since the results obtained correspond to noninteracting 
particles they can be also explained in a single particle 
picture due to the realization of Bragg conditions, and are 
also valid for bosons. We have studied the consequences 
of the previous confinement in the nonequilibrium dy- 
namics of trapped particles in ID when the center of the 
trap is suddenly displaced, and confirmed evolution of 
the center of mass obtained in recent experiments. 

The region over which particles are confined in the trap 
can be controlled in various ways. The most obvious one 
is by changing the strength of the confining potential, 
where the extension of such regions can be regulated. 
Other way is changing the periodicity of the lattice, which 
leads to a different "slicing" of the system. The change 
of the periodicity also generates in the fermionic case 
the possibility of obtaining local insulating phases with 
sizes that can be controlled changing the strength of the 
additional alternating potential. This gives rise to a pic- 
ture that is similar in some aspects to the the Hubbard 
model analyzed in Rcfs. [23, El- We have shown that 
although insulating phases appear in this nonintcracting 
case, the gaps that are locally opened are not seen in the 
single particle spectrum. In order to observe them it is 
necessary to study the local density of states. The local 
compressibility defined in Refs. |20t |21| was also proven 
to be a genuine local order parameter to characterize the 
new insulating phases since it is always zero there. A 
scalable phase diagram for these systems was also pre- 
sented. Finally, we considered the two-dimensional case 
and the formation of insulating regions due to the pres- 
ence of periodic potentials. We showed that the local 
compressibility also characterizes those 2D regions in an 
unambiguous way. 
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